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OPTIMIZING INTERMITTENT WATER SUPPLY IN URBAN PIPE 
DISTRIBUTION NETWORKS* 

ANNA M. LIEBtt, CHRIS H. RYCROFTS*, AND JON WILKENINGtt 


Abstract. In many urban areas of the developing world, piped water is supplied only intermit¬ 
tently, as valves direct water to different parts of the water distribution system at different times. 
The flow is transient, and may transition between free-surface and pressurized, resulting in complex 
dynamical features with important consequences for water suppliers and users. Here, we develop a 
computational model of transition, transient pipe flow in a network, accounting for a wide variety 
of realistic boundary conditions. We validate the model against several published data sets, and 
demonstrate its use on a real pipe network. The model is extended to consider several optimization 
problems motivated by realistic scenarios. We demonstrate how to infer water flow in a small pipe 
network from a single pressure sensor, and show how to control water inflow to minimize damaging 
pressure transients. 


1. Introduction. From the dry taps of Mumbai to the dusty reservoirs of Sao 
Paolo, urban water scarcity is a common condition of the present, and a likely feature 
of the future. Hundreds of millions of people worldwide are connected to water 
distribution systems subject to intermittency. This intermittent water supply may 
take many forms, from unexpected disruptions to planned supply cycles where pipes 
are filled and emptied regularly to shift water between different parts of the network 
at different times [17, 33]. In Mumbai, for example, Vaivaramoorthy [33] reports 
that on average, residents have water flowing from their taps less than 8 out of 
24 hours. Intermittent supply is often inequitable, with low-income neighborhoods 
experiencing lower water pressure and shorter supply durations than high-income 
ones [31]. Intermittent supply not only limits water availability, but also compromises 
water quality and damages infrastructure. With field data from urban India, Kumpel 
and Nelson [18] quantified the deleterious effect of intermittency on water quality, 
showing that both the initial flushing of water through empty pipes—as well as 
periods of low pressure—corresponded with periods of increased turbidity and bacterial 
contamination. Christodoulou [7] observed when that a drought in Cypress ushered in 
two years of intermittent supply, pipe ruptures increased by 30%-70% per year. 

Whereas intermittent water supply creates challenges for water managers and 
water users, the phenomenon creates opportunities for applied mathematics. It is an 
interesting and difficult mathematical problem to efficiently model transient pipe flow 
in networks—including transitions to and from pressurized states—with uncertain or 
complex boundary conditions. In this work we introduce a framework to not only 
describe intermittent water supply, but also use optimization to improve either our 
description of the system, or the operation of the described system in order to reduce 
risks such as infrastructure damage. 

Intermittent supply falls in somewhat of a modeling gap. Water distribution 
software abounds, including the free and open source software EPANET [27] produced 
by the US government, as well as many commercial packages [13]. Yet, to the 
authors’ knowledge, all these fail to account for filling, emptying, and instances 
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of subatmospheric pressure—phenomena that are vitally important for users and 
managers dealing with intermittent supply. Sewer system software such as the Storm 
Water Management Model (SWMM) [28] and Illinois Transient Model (ITM) [20] to 
some extent handle the physics of interest, but are packaged in elaborate graphical 
user interfaces and are not readily amenable to model improvements or optimization. 

Furthermore, the authors have encountered a relative paucity of research work 
dealing with modeling intermittent supply. The work of De Marchis [9] explicitly 
studies filling and emptying in a water distribution system in Palermo, Italy, but with 
a method of characteristics implementation of the classical water hammer equations. 
This treatment assumes pipes are either entirely dry or entirely full, and that air 
pressure inside the pipes is always atmospheric. After calibrating a friction parameter, 
they found about 5% agreement with empirical data. Subsequent work reported by De 
Marchis [10] uses the same model to assess losses in the distribution system. Freni [12], 
uses this model to determine pressure valve settings to reduce distribution inequality, 
but through scenario comparison rather than optimization. For sewer flow, Sanders [30] 
presents a network implementation of the two-component pressure approach (TPA) of 
Vasconcelos [37]. The modeling for ITM was published by Leon in [23]. Urban water 
drainage is coupled to free surface flow by Porsche and Klaar [1]. Note that Buosso et 
al. [5] give a general review of the storm water drainage literature with more details 
than we have provided here. 

The present work comprises an effort to address the scarcity of tools available for 
those interested in modeling the details of intermittent supply, and to specifically incor¬ 
porate such tools within an optimization framework. We use an underlying model of 
coupled systems of one-dimensional hyperbolic conservation laws that strikes a balance 
between real-world relevance and both computational and theoretical tractability. Our 
computational framework will allow for straightforward implementation of alternative 
physical models in future studies. 

2. Model. The Preissman slot formulation [25] is used to describe flow within 
each pipe, building on existing literature for transient, transition flow in closed conduits. 
The flow is assumed to be inviscid and incompressible. The dynamical description 
considers depth-averaged flow within a modified geometry that permits a single set 
of equations to describe both free-surface and pressurized flow. Consideration of 
one-dimensional dynamics is a reasonable approximation given that the ratio of pipe 
diameter D to pipe length L is 1% or smaller in realistic scenarios. 

The modeled variables are the cross-sectional area A and the discharge Q, which 
can be used to compute quantities of practical interest such as pressure head and 
velocity. Pressurization effects, though features of compressible flow, are accounted for 
by the Preissman slot cross section pictured in Figure 2.1(a). Above a transition height 
Uf near the physical pipe crown, water fills the narrow slot, contributing hydrostatic 
pressure that may be interpreted as the pressure within the full pipe. The slot width 
is related to the effective pressure wave speed a in the pressurized pipe via = gAf /Tg, 
where g is the acceleration due to gravity and A/ is the cross-sectional area at the 
transition height yf (determined by the pipe diameter and Tg ). The pressure wave 
speed a is, to first order, a function of pipe material only. In practice, the slot width 
Tg is a parameter chosen to set the value of pressure wave speed a, which may range 
from 20-1250 m/s [37] depending on practical context and numerical constraints. In 
each pipe, the governing equations for A and Q are the de St. Venant equations for 
free-surface flow [16, 23], 

qt -I- (F(q))a; = S, 0 < X < L, 0 < t < T, 
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Fig. 2.1. (a) The Preissman slot model describes flow in a pipe of diameter D by depth-averaging 
over this cross-section. Each filled cross-sectional area A corresponds to a height ^ and length l{i). 
When ^ exceeds a transition height yf, water fills a narrow slot with width Ts and contributes 
additional hydrostatic pressure, (b) An example pipe network. 


where L is the length of the pipe, T is the duration of the study period, and 


q = 



F = 


Q \ 



( 2 . 2 ) 


where g is the acceleration due to gravity and the pressure contribution I {A) is given 

by 

rh{A) 

I{A)= {h{A)-Omd^, (2.3) 

Jo 

where l{^) is the pipe width at height The quantity p = gI{A)/{A) is the average 
hydrostatic pressure over a cross-section. In the Preissman slot description, the 
pressure head H is entirely captured by the hydrostatic term p, and given by 


~ P9 A 

The friction term S is 


(2.4) 


S = (So- Sf)gA, 


(2.5) 


where So is the slope of pipe bottom and S/ empirically accounts for friction losses. 
In what follows we use the Manning equation 


^ A^Rh(A)'^/^ 


( 2 . 6 ) 


where Rh(A) = AjP^ is the hydraulic radius, which depends on the wetted perimeter 
P-u], which is a function of A and D. The constant Mr is the Manning roughness 
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coefficient, which has an empirical value depending on pipe material. Other empir¬ 
ical formulae such as Hazen-Willams or Darcy-Weisbach may also be used. Initial 
conditions q(a;, 0) are assumed to be known, and boundary conditions q(0, t) and 
q(L, t) are assigned based on the network connectivity and external inputs described 
in Section 3.2. 

The Preissman slot approach has been used, for example, by Trajkovic et al. [32], 
who compared the model with experimental data; by Kerger et al. [16], who imple¬ 
mented an exact Riemann solver and introduced a “negative slot” modihcation to 
handle subatmospheric pressure; by Leon et al., [22], who allowed the pressure wave 
speed to vary slightly; and by Borsche and Klar [1] , who used a hexagonal cross-section 
to simplify the area-height relationship. Other transient flow models for air and water 
in single pipes include a “two-component pressure” approach [37, 34]; a single-equation 
model with a modified pressure term [3, 4]; a two-component model [23]; and a 
three-phase model accounting for air, air-water mixture, and water [15]. 

Note that the de St. Venant equations themselves make no assumptions about the 
channel cross section. The numerical methods in Section 3 may be rather easily adapted 
to other cross-sectional geometries, and in our implementation we also include the 
option to simulate flow in channels with uniform cross-sections. Our implementation 
could therefore be used to simulate other networks like rivers or irrigation canals, but 
such options have not yet been explored. 

3. Numerical Methods. 

3.1. Pipes. Each pipe k is assigned a local coordinate 0 < Xk < Lk and divided 
into Nk cells of width Axk, where the cell centers are at Xkj = j + \Axk and the cell 
boundaries are at Xk^j±i /2 = Xk,j ± \ Axk for j = 0,..., Nk — 1, as shown in Figure 
3.1. Each pipe is also padded with ghost cells at Xk,-i and Xk,Nk that are used to set 
the boundary conditions, described in Section 3.2. A total of M timesteps of size At 
are taken. 



Fig. 3.1. Spatial grid layout for finite volume method in pipe k of length Lf^. 

In what follows we drop the k subscript and assume we are working within a single, 
specified pipe. Cell averages of area and discharge (A", Q") = q" at the timestep 
are updated with an explicit third-order Runge-Kutta total variation diminishing 
(TVD) scheme [14] that may be written in terms of first-order Euler update steps £’(q) 
as 

q=^q-+^-E{E{qn), q"+' = ^q" + ^^^(q). (3.1) 
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Each Euler step E{q) consists of updating the conservation law terms, the source 
terms, and the ghost cell values. For the interior cells, the update is 

E(q)=E,(E,(q)) (3.2) 

where the subscript c denotes the conservation law update and the subscript s denotes 
the source term update. First we treat the homogeneous conservation law term 

q* + (F(q)), = 0 (3.3) 


with a Godunov update 

For the numerical flux function F, we use a Harten-Lax-van Leer (HLL) Riemann 
solver similar to that of Leon et al. [22], which approximates the solution to the 
Riemann problem between a pair of cells (with left and right states q^ and q^j 
respectively) with a center state q* separated from the left and right states by shocks 
with speeds sl and sr, respectively [24]. 

The expression for the center state flux F* = F(q*) is found by applying the 
Rankine-Hugoniot condition across each shock. The Godunov update for this scheme 
is found by sampling this solution structure to obtain 


where 


Fl = F(qL) if Sl > 0, 

F = if Si < 0 < Sr, 

Ffl = F(qi;) if Sij < 0, 


(qL,qfl) 


(qi,qi+i) atj + i, 
(qj-i,qj) atj-i. 


We computed the shock speeds via 


(3.5) 


(3.6) 


Sl = UL — ^L, Sr = Ur + ^Ir 


where Uj = Aj/Qj is the velocity in cell j and 


( l igI(A.)-gIiAAY^ 

n,- = J V A,{A,-Aj) 

lc(A*) 


if A,,i > Aj T e, 
if A^ ^ Aj “t“ c, 


(3.7) 


where c{A) = sfgAfli^ is the gravity wave speed [22] . The center state is approxi¬ 
mated by linearizing the equations to obtain 




(3.8) 


where c = {ul + ur)/2. Note that if Aj < A*, the separating wave is not in fact a 
shock, but a rarefaction, and the expression for Sj gives the speed of the head (or tail) 
of the appropriate rarefaction wave [21]. A tolerance of e = 10“® is incorporated into 
(3.7) to account for the possibility that when A* and Aj are very close, small errors in 
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the evaluation of I {A) may cause the expression underneath the radical to become 
negative. By Taylor expanding the expression underneath the radical at Aj^ one can 
verify that the two cases in (3.7) connect continuously, and thus the incorporation of e 
has a minimal effect on the computation. The shock speeds sl and sr may also be 
estimated by evaluating u and c{A) for Roe averages A and u and taking minima and 
maxima for left and right waves, respectively [30], but the authors found the dynamics 
of interest in the present work were less robust with this treatment. 

Computing the HLL fluxes requires frequently evaluating the pressure term I {A), 
the wave speed c{A), and the integral which are not analytic expressions of A 

for the Preissman slot geometry. To avoid rootfinding at every cell at every time step, 
Chebyshev polynomials were used to accurately express these functions of A and their 
inverses. Several of the functions are ill-suited to polynomial approximation due to 
fractional power singularities at either end of their domain, so we implemented an 
accurate interpolation by expanding in a series involving the relevant fractional power 
(Appendix A). 

The friction and slope source term updates, which only affect the second component, 
take the form 

Es{Q^) = Q” + At 5 -b • (3.9) 


3.2. Junctions. Each pipe domain is padded with ghost cells, which serve to 
compute fluxes to update boundary cells in each pipe’s computational domain. In what 
follows, we use the notation (A^^j, to denote the ghost cell values of A and Q at 
the relevant end of pipe k and denote the values of the last cell in the computational 
domain by The algorithm uses the network layout to determine what 

update routine is used on the ghost cell values (Ag,^^, The update routines 

fall into three categories: external boundaries, two-pipe junctions, and three-pipe 
junctions. Note that the user must specify additional information for the external 
boundary routine. For other boundary treatments in networks, see [30], [6], and [19]. 

3.2.1. External boundary routine. When one end of a pipe is connected to 
the outside world, the user may specify a variety of cases, summarized in Figure 3.2, 
to describe the external conditions. In each case, the user specification of boundary 
condition type allows the solver to update the ghost cell values. For example, in the 
network shown in Figure 2.1(b), a time series of Q or A would be specified at each 
“inflow” node (cases (2) and (3)). At the labeled “outflow” nodes, the user may choose 
between several possible descriptions of how a valve allows water to exit the end of a 
pipe. These boundaries may be considered either as either orifice outflow (case (4)) or 
unimpeded outflow where no waves are reflected (case (1)). A closed valve may be 
simulated as a reflective boundary (case (0)). 

As only one pipe is involved, in what follows we drop the k superscript on the 
ghost and internal cell values. For cases (0) and (1), the ghost cell values (Aext,Qext) 
are updated via 


( Aext; 


Qext) — 


(^ in; Qin) 
(^ in; Qin) 


case (0): reflect everything, 
case (1): reflect nothing. 


(3.10) 


Case (0) reflects all waves and case (1) reflects none [24]. Physically, case (0) cor¬ 
responds to a dead end and case (1) corresponds to an opening with unimpeded 
outflow. 
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Fig. 3.2. Single pipe boundary cases. User may specify either reflection of all or no waves, a 
value for one of the dynamical variables A or Q, or a valve opening. 


Cases (2) and (3) arise when the user specifies a time series for either Aext{x,t) or 
QextiXyt) at the external boundary. During each Euler update step, the solver first 
determines whether the interior flow is super- or subcritical, depending on whether 
the Froude number Fr = u/c{A) is greater than or less than unity, respectively. 

In the supercritical case (2), if Qext < 0 at x = 0 or Qext > 0 at x = L, then 
outflow case (2.1) applies. Information from the boundary cannot propagate inside the 
domain under these conditions, and thus case (2.2) is evaluated in the same manner 
as the extrapolation case (1), where all outgoing waves continue untrammeled. For 
supercritical inflow, case (2.1) applies. If Qext is specified, the undetermined ghost 
cell value of ^ext is set to the initial inflow cross-sectional area; otherwise the Froude 
number in the ghost cell is set equal to the Froude number just inside the domain. 

For subcritical flow, information may propagate in both directions, and the solver 
determines a value for the unspecified component ^ext or Qext- The approach in the 
current work is to attempt to follow an outgoing characteristic and use the value of the 
Riemann invariant along this characteristic to solve for the unknown external value. 
Sanders and Katopodes [29] used an exact version of this for networks of channels with 
uniform cross-section (in this case the Riemann invariants are simple). The same idea 
was implemented iteratively for a circular geometry by Leon et al. [21] and Vasconcelos 
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et al. [37]. However, care must be taken to ensure that the characteristic assumption 
is valid. Thus our algorithm for case (3.1) works as follows: 

1. Calculate the Riemann invariant for the last cell in computational domain 
according to = QinMin ± where (Ain, Qin) are values in the last cell, and 

4> = J{c/A)dA depends on geometry. (j){A) = 2^^gA/w for uniform cross-sections of 
width rc, and has no analytic expression for circular cross sections (see Appendix A 
for Chebyshev representation). 

2a. If Aext is specified, set 

Qext — Anxt (Qin/Ain ± WAn) “ '('(A ext ))), 


where the solution is physically valid since Q is allowed to be positive or negative. 
Even though a solution will always be found, it may violate the subcritical condition. 

2b. If Qext is specified, check for compatibility as described below. If incompatible, 
go to case (3.2). Otherwise, rootfind to solve for Aext satisfying 

Qin/Ain ± /'(Ain) = Qext/Agxt i /)(Aext). 


3. Update ghost cell values to (Aext,Qext)- 
Rootfinding in step (2b) above is by no means guaranteed to work; indeed, a positive 
solution Aext only exists for certain combinations of the parameters. Physically, this 
means that not all values of Qext may be continuously connected to the interior 
state, and choosing certain Qext forces a shock between the ghost cell and the last 
computational cell. For the boundary at x = 0, for a given value of Qext < 0 
there is a maximum allowed value of c_ = Qin/Ain — (/(Ain). Similarly, for the 
boundary at x = L, for a given Qext > 0, there is a minimum allowed value of 
c+ = Qin/Ain -I- (/(Ain). For the uniform cross-section case, one can find analytic 
formulas for these maximum/minimum allowed values: 




(3.11) 


For the Preissman slot, to find the values of Cj- we estimate the critical point of the 
function g(t) = Qext/^ ± (/(i) rather than rootfind for the exact value. We define 
Cj- = Qext/® ± </(®), where = yQext pip® diameter D), which corresponds to 
approximating the gravity wave speed by y/gh{A). The compatibility condition is 
then 


if 


f c_ > c* 1 

1 ] 

\ x = 0 1 

{ 

r { 


\ c+ < c; j 

1 1 

1 - = ^ i 


then Qext is not compatible. 


(3.12) 


Should Qext fail the compatibility condition, case (3.2) applies, and the solver 
sets Aext = Ain. This choice causes the HLL center state estimate A* < Ajn = Aext, 
consistent with states separated by shocks in the eyes of the approximate Riemann 
solver update. 

Lastly, if a valve or gate opening is specified, case (4) applies. Bernoulli’s equation 
applied to flow through an orifice gives 

Qext — CdA\/2g{h\n Cch ext ), (3.13) 


where hin = /i(Ain), /text is a height between 0 and D that indicates the valve opening, 
g is the acceleration due to gravity, and A = A(/iext) is the area of the outflow orifice. 
The discharge coefficient Cd = 0.78 and the contraction coefficient Cc = 0.83, are 
empirical constants from experiment [32]. In this case, Aext = Ain. 






Fig. 3.3. Junction routine schematics for (a) two pipes, and (b) three pipes. 


3.2.2. Two-pipe junction routine. In the absence of valves, we apply mass 

conservation and assert that water height is constant across a two-pipe junction. Hence, 
the ghost cell in pipe ki is updated by translating water height from the pipe ^2 to 
the local geometry in pipe ki and copying the discharge from pipe /c 2 . That is, for 
example, pipe fci gets ghost cell values (^ext) Qeit) such that = hk^{A\l) 

and Qext = Qfn- This procedure is shown in Figure 3.3(a). 

3.2.3. Three-pipe junction routine. Triple junctions are divided into three 

pairs of two-pipe subproblems, which are solved to find fluxes. Each two-pipe sub¬ 
problem contributes half the flux to an incoming pipe, as shown in Figure 3.3(b). 
For example, the fluxes at the end of pipe ki are the sum of half the flux from the 
two-pipe junction routine between (ylfj(,(5fn) and and half the flux from 

the two-pipe junction routine between (^fn’Qfn) coordinate 

systems do not point in the same direction (e.g. both pipes have a; = 0 at the triple 
junction), then one of the discharge terms must have a relative minus sign before the 
two-pipe junction routine is applied. The flux assignment does not otherwise depend 
on geometry. 

The authors believe this formulation is a simple way to couple the one-dimensional 
problems with minimal computational effort, since it recycles the two-junction solver 
and introduces no new data structures or solution routines. This routine may be 
improved by accounting for geometric effects and energy losses (often referred to as 
“minor losses” in pipe flow parlance) due to the specific geometry of the junction. Other 
treatments of these types of boundaries include Leon et al. [19], in which the triple 
junction has hnite area and a dropshaft, and n-pipe junction implementations [8, 1]. 

3.3. Model implementation. The computational model (provided in supple¬ 
mentary materials) is written in C-|—h, and the different simulations are initialized 
through two text files: an EPANET-compatible file describing network layout, and a 
configuration file specifying additional simulation parameters. In addition, a Cython 
wrapper is available, allowing simulations to be launched from iPython notebooks. 
The simulation running time analyses that are reported in the following sections were 
performed on a Mac Pro (Mid 2010) with dual 2.4 GHz quad-core Intel Xeon processors. 
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Fig. 4.1. Model (thin lines) and experimental [32] (points) pressure head at locations P7 and 
P5, for gate openings eo = 0.008 m (top) and eo = 0.015 m (bottom). 


using a single thread unless stated otherwise. The authors found that parallelizing 
the network using the OpenMP library (so that all pipes are solved simultaneously) 
offered no advantage, because doing single-timestep Riemann solver updates along 
each pipe is too fast to make it worthwhile to spawn separate threads. 

4. Model validation and results. 

4.1. Experiments of Trajkovic et al. (1999). Experiments were performed in 
a single pipe of flow transitioning from open-channel to pressurized and back [32]. The 
experimenters used the Preismann slot to model their results, and their experiments 
have been revisited in several other studies [22, 36]. The experiments were performed 
in a single pipe of length L = 10 m and diameter D = 0.1 m, set at a slope of 
2.7% as reported in [32]. The Manning roughness coefficient was estimated as 0.002 
after examining modeling results for a small range of roughness coefficient values. 
We consider the “type A” experiments where the initial condition was supercritical 
unpressurized flow with h = 0.014 m and Q = 0.0013 m/s throughout the pipe. At 
time t = 0 s a gate at a: = 9.6 m is closed, causing the pipe to pressurize. At time 
t = 30 s, the gate is partially reopened to a height eo, which allows the flow to de¬ 
pressurize if eo is sufficiently large. This experiment was modeled with Ax = 0.05 m 
and a Courant number of 0.6. Pressure traces at two sensors, P5 and P7, located at 
0.5 m and 2.5 m upstream of the gate, respectively, are plotted in Figure 4.1 for three 
values of eg are shown below. Comparison with experimental data shows that the 
model arrival times and pressure traces are in good agreement. 

4.2. Water Hammer. In this textbook example, similar to a case presented by 
Vasconcelos et al. [36], a sudden valve closure gives rise to a so-called water hammer 
phenomenon in a pressurized pipe, whereby a shock wave of high pressure travels 
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Fig. 4.2. Water hammer example. Time series for H{x = L,t)/Hg and Q{x = 0,t)/Qg. 


rapidly up and down the pipe. The setup for this problem consists of a horizontal 
frictionless pipe of length L = 600 m and diameter D = 0.5 m. Initially, the pipe has 
pressure head Hq = 150 m, and a discharge Qq = 0.1 m^/s. The upstream pressure 
head at a: = 0 m is held steady at 150 m, and a reflective boundary is applied ai x = L. 
The pressure wave speed set by the Preissman slot is 1200 m/s. Results in Figure 4.2 
show the time evolution of the nondimensionalized pressure head H{x = L,t)/Ho at 
the valve and the discharge Q{x = Q,t)/Qo at the upstream end. The modeled results 
agree with the classical water hammer equations [40], which assert that for a pressure 
wave propagating in pressurized pipe, the relationship between the change in pressure 
head AH and the change in velocity AV = AQjAj is 

AH _ a 

where Af is the cross-sectional area of the full pipe (corresponding to transition height 
Uf), g is the acceleration due to gravity and a is the pressure wave speed. With 600 
grid cells and a Courant number of 0.6, we computed \ ~ gl ~ 0.00002. 

4.3. Grid refinement. We consider two horizontal, frictionless pipes connected 
in serial, each of length 500 m. The left pipe (denoted pipe 0) has diameter 1 m, 
and the right pipe (denoted pipe 1) has diameter 0.8 m, allowing for examination of 
the effect of a constriction upon the simulated flow. Starting with initial conditions 
Q{x,t = 0) = 1 m^/s and h(x,t = 0) = 0.75 m in both pipes, we run the simulation 
until t = 100 s. The boundary conditions are specified inflow Q{x = 0,t) = 0.1 m^/s 
for pipe 0 and and reflection at x = 500 m for pipe 1. Pressurization occurs in part of 
the system, demonstrating the Preissman slot in action. Note that these examples 
used a pressure wave speed of 12 m/s corresponding to a relatively wide slot. Running 
convergence studies with a physical wave speed on the order of 1000 m/s would have 
given the same behavior but requires a much smaller timestep. We run the simulation 
with a range of values of Ax to obtain the results pictured in Figure 4.3, which shows 
final pressure head h{x,t = 100) in both pipes. Table 4.1 shows the error at time T, 
defined as ||e||i = k 1^^ ~where h^{xj,T) is the solution in pipe k on 
the finest grid, evaluated at point Xj. As expected, we observe Hrst-order convergence, 
due to the presence of shocks. The table also lists the computation times of the test. 
The computation times scale quadratically with N since the number of timepoints M 
must also scale with N to satisfy the CFL condition. 

4.4. Realistic network. We next consider, as an illustration of the model’s 
ability to handle a realistic scenario, a larger network layout based on a network 
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Fig. 4.3. Grid convergence test, water height h(x, t) at t = 100 s. 
Table 4.1 

Convergence results, solution computed at time t = 100 s. 


N 

M 

Wall clock time (s) 

l|e||l 

50 

200 

0.10 

0.18251 

100 

400 

0.40 

0.10791 

200 

800 

1.54 

0.05608 

400 

1600 

5.84 

0.02188 

800 

3200 

23.80 

0.01380 


experiencing intermittent supply in a suburb of Panama City. The network contains 
fourteen pipes, ranging 300-1100 m in length and 0.4-1.0 m in diameter. Slope terms 
are non-zero in all pipes and Manning coefficient is estimated at 0.015 based on pipe 
material. We simulate filling over a period of 1120 seconds, starting with initial 
conditions pictured in Figure 4.4(a), to obtain the pressures in Figure 4.4(b)-(f). The 
simulated pressure values fall within in a realistic range, as do the filling times of the 
pipes that pressurize. The fact that not all pipes pressurize during the simulation 
period is also commensurate with observations. 

5. Optimization. We now extend the computational model to consider several 
optimization problems that are inspired by realistic scenarios. We use the Levenberg- 
Marquardt (LM) algorithm, a trust-region method, to determine an unknown quantity 
that minimizes an objective function. When the unknown quantity is a time series, 
as in Sections 5.1 and 5.2, the degrees of freedom describe either a Hermite spline or 
a truncated Fourier series. The LM algorithm requires computing a Jacobian, and 
this was calculated using finite differences, since the nonlinearly coupled boundary 
conditions made a variational formulation impractical. The Jacobian columns are 
computed in parallel using the OpenMP library, which provided considerable time 
savings over a serial implementation for the problems considered. 

5.1. Recover unknown boundary data. The uncertainty of intermittent wa¬ 
ter supply motivates our first example, in which we are interested in recovering 
unknown boundary information based on our knowledge of pressure somewhere in the 
network. Below, we provide proof-of-concept for the modeling/optimization framework 
by solving a problem where we know the correct answer. Figure 4.5 shows a diagram 
of the example network considered. The pipe diameters are all 1 m, and the lengths 
are 500 m, 500 m, and 125 m for pipes 0, 1, and 2, respectively, Aa; = 2.5 m for 
all pipes, and the number of time steps is M = 700. Inflow to pipe 0 is a constant 
1 m^/s. Pipe 1 has a reflective boundary at the downstream end, and the outflow time 
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Fig. 4.4. Pressure head (m) in example (4-4)■ Grey pipes are empty. Times shown are (a) 
t = 0 s (h) t = 2A0 s, (c) t = 480 s, (d) t = 720 s, (e) t = 960 s, (f) t = 1200 s. 


Reflect 



Fig. 4.5. Illustration of boundary recovery example. Inflow Q(0,t) to pipe 0 is prescribed on the 
left. Measurements of h{x* fl) = h*{t) in pipe 0 are used to find the unknown pipe 2 outflow Q{L,t). 


series of pipe 2 is not known. We generate a test data set with the ontflow time series 
QtTVie{L,t) as depicted in Figure 5.1(a). The slot width is 0.00053 m, corresponding 
to a pressure wave speed of 120 m/s. where H{x*,ti) is the solution at the sensor at 
time step i. The degrees of freedom are Hermite spline coefficients or Fourier modes 
for the boundary value time series Q{L,t) on the outflow end of pipe 2. 

We initialized the optimization routine with Qi{L,t) = 0. The time series H*{t) 
at the sensor resulting from the true boundary condition Qtrue and the initial proposed 
boundary condition Qi are shown in Figure 5.1(b). The optimized time series Qf 
found with Hermite and Fourier representations are pictured against the initializing 
and true time series in Figure 5.2(a). The run times, ratio of initial objective function 
value fi to final value /j, and error in Q(L,t) before and after optimization are shown 
in Table 5.1. We considered both Hermite splines and Fourier modes with 8 degrees 
of freedom. Results are shown in Figure 5.2. Notice that after time ST = a/5x, 
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Fig. 5.1. (a) Pipe 2 outflow times series Qtrue used to generate data, (b) Sensor measurement 
H*{t) when pipe 2 outflow is given by Qtrue (black dashed line) vs. Qi (red dash-dot line). 




(a) 


(b) 


Fig. 5.2. (a) Recovered time series Qf compared with Qtrue and Qt. (b) Error 

\\Q{x,ti) — Qtrue{x,ti)\\ at time tt, where || ■ ||2 is the discrete L 2 norm over the spatial domain. 


information has not had time to propagate from the boundary to the sensor so the 
Hermite representation has difficulty thereafter. The Fourier series representation 
maintains a low level of error through the whole time series as the true solution has 
the same value at t = 0 as at t = T. 

5.2. Reduction of variation in pressure. We next consider controlling bound¬ 
ary conditions to reduce potential damage to pipes in the water distribution system. 
Although the causes of pipe damage are complex and myriad, pressure transients play 
a role [11, 2]. Internal pressure contributes to axial, longitudinal and hoop stress 
in pipes [26], and pressure transients may help cause both pipe bursts (associated 
with high-pressure events) and collapses (associated low-pressure events) [39]. During 
filling, the presence of trapped air may amplify pressure peaks and further increase 
damage [35, 18]. Within the scope of the Preissman slot model, which does not 
explicitly describe air entrainment, we consider the potential for pipe damage by 
measuring, at each time step ti, an estimated total variation in the pressure head H, 
which we denote {dH/dx)i, defined as 

K Nk-l 

{dH/dx), (5.1) 

k—1 j—0 

where K denotes the number of pipes in the network, Np the number of grid cells in 
pipe k, M the number of time slices, and Wj. - = H{A{xj,ti)) for A in pipe k. Note 
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Table 5.1 

Results for boundary data recovery. 



Hermite 

Fourier 

CPU time (s) 

176 

135 

wall clock time (s) 

64 

51 

parallel speedup 

2.8 

2.6 

h 

6.9 X lO'^ 

6.9 X 10"“ 

ff 

1.9 X IQ-i 

3.2 X 10-2 

h/ff 

2.7 X 10“® 

4.6 X 10“'^ 

\\Qi Qtrue 11 

22.6 

22.6 

11Q/ Qtrue 11 

0.477 

0.296 


Table 5.2 

Mean and maximum {dH/dx) for different L 1 /L 2 . 


L 1 /L 2 

Mean (dH/dx) 

Max (dH/dx) 

0.25 

5.0609 

48.3832 

0.50 

5.4465 

60.9356 

0.75 

6.1335 

60.6128 

1.00 

3.4518 

24.4900 

1.25 

5.5090 

40.9716 

1.50 

5.1254 

34.4893 

1.75 

5.1019 

34.5001 

2.00 

5.1504 

34.6804 


that 


(dH/dx) 



dH 

lo 

dx 


dx 


(5.2) 


where dH/dx may only exist in a weak sense due to the presence of shocks. The 
quantity (dH/dx) serves as proxy for pressure gradient that penalizes large water 
hammers associated with both high- and low- pressure events, and its minimization 
allows for exploration of operational regimes that potentially alleviate hydraulic 
contribution to pipe damage. 

As we are interested in studying existing networks, rather than planning new ones, 
we first study how network geometry affects the time evolution of (dH/dx). We use 
the network structure shown in Figure 4.5 to study how the relative lengths of pipes in 
the network affects the interaction of reflected waves coming through a triple junction. 
We fix Lq = Li = 100 m, and L 2 is varied between 25 m and 200 m. All boundaries 
have (5 = 0, and the initial conditions are h = 0.8 m in all pipes and Q = 2 m^/s in 
pipe 0, and Q = 1 m^/s in pipes 1 and 2. These conditions were chosen to give a 
scenario where the flow transitions from free-surface to pressurized. The simulation 
time is 18 s and the pressure wave speed is 120 m/s. A time series of (dH/dx) for 
each network geometry is shown in Figure 5.3. 

Mean and maximum values are tabulated in Table 5.2. Interestingly, for Li = L 2 
we see the least dramatic pressure variation. Although the peak value of (dH/dx) 
varies by up to a factor of two over these different geometries, the mean values are 
quite comparable, suggesting that boundary control may still be useful even in different 
network geometries where reflected waves interact differently. 

We now use the optimization framework to suggest inflow patterns minimizing 
(dH/dx) for one of these networks under a pressurization scenario. We define the 
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Fig. 5.3. Time series of {dH/dx) for different LifL 2 . 


objective function / = i||r|p, where the components of r are given by 


n = 



{dH/dx)i 


(5.3) 


for i = 1,..., M. Motivated by practical applications, we suppose we want to supply 
a fixed volume of water Vin over a time period T, and vary the inflow time series 
(5(0, t) at the left end of pipe 0 in order to minimize /. Below we present three results 
highlighting both the usefulness and limitations of the framework. In all cases, the 
network connectivity is kept the same, and boundaries for the outflow ends of pipes 
1 and 2 are set to (5 = 0. The initial conditions, time series representation, and 
simulation time are varied as follows: 

(I) Lq = Li = 100 m, and L 2 = 25 m, roughness coefficient = 0.015. Initial 
condition is all pipes empty. Degrees of freedom are cubic Hermite spline coefficients, 
with the first point Q{0, 0) determined by setting the integral of the Hermite time 
series equal to the desired inflow Vjn. The initial constant inflow is compared with the 
optimized time series in Figure 5.4(a). How {dH/dx) varies in time is shown in Figure 
5.4(b). 

(H) Lq = Li = 100 m, and L 2 = 50 m, roughness coefficient = 0.008. Initial 
condition is Q{x,t) = 5 m/s in pipe 0, with pipes 1 and 2 empty. The degrees of 
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Table 5.3 

Numerical results for pressure gradient reduction. 


Case 

DOF 

a 

T 

M 

CPU time 

Wall time 

fo 

ff/h 

Fin 

I 

15 (Hermite) 

100 

60 

15000 

5298 

772 

13.3013 

0.2211 

120 

II 

15 (Fourier) 

10 

120 

4000 

2819 

405 

1.53498 

0.7401 

150 

III 

15 (Fourier) 

100 

16 

3600 

2576 

350 

1.83238 

0.9988 

20 


freedom are Fourier modes, with the constant mode constrained to 2Vi^lT, such that 
Vin is the total inflow volume. The optimization starts with a time series Q{0,t) for 
pipe 0 which approximates a step function inflow. Pressure wave speed a = 10 m/s 
and simulation time T = 120 s. Results are pictured in Figure 5.5. 

(Ill) As in case (II), but with pressure wave speed a = 100 m/s and T = 16 s. 
Results for case (II) were computed with a relatively wide slot width, corresponding 
to a gravity wave speed of a = 10 m/s. With a more realistic pressure wave speed 
of 100 m/s in case (III), the algorithm found only very slight improvements for the 
inflow time series, as shown in the table above. This is probably due to the presence 
of unphysical oscillations accompanying the initial pressurization event. That the 
Preissman slot model, with a narrow slot width, gives rise to such oscillations during 
rapid pressurization has been noted previously by Vasconcelos [38], among others. 
Vasconcelos suggests filtering or artificial viscosity to damp out these fluctuations, 
but also concedes the averaged behavior may present a reasonable approximation of 
reality. Even if these oscillations are filtered, averaged, or damped into submission, 
they limit the utility of the Preissman slot model for informing smooth optimization 
when rapid pressure transients are present. However, in regimes with more gradual 
filling or emptying, the Preissman slot model may suffice for both pressure gradient 
optimization and boundary value recovery. 

6. Conclusions and Future Work. In this work we have introduced the study 
and management of intermittent water supply as a problem in applied mathematics, a 
problem with interesting theoretical facets as well as compelling practical implications. 
We present results from a combined modeling and optimization framework that captures 
some of the important dynamics in the systems of interest and demonstrates two 
optimization examples with real-world significance. Although the present results 
may help guide some ideas for managing pressurized flow, model improvements and 
alternative optimization schemes may be necessary to address applications of pressure 
gradient reduction and other objectives in larger and/or poorly-mapped networks. 

In the future we hope to more carefully extract the physics and timescales relevant 
to the application of intermittent water supply in a robust model that retains the 
computational tractability of the current one. The code allows for straightforward 
implementation of alternative models for single-pipe flow and boundary coupling, 
allowing us to make the pressure gradient reduction optimizations more relevant to 
real-world data. We also hope deal with questions of water quality and to extend our 
implementation to tackle larger networks and other objective functions, aspiring to 
examine, understand, and perhaps improve real-world scenarios. 

Acknowledgments. We thank John Erickson and Kara Nelson (UC Berkeley 
Civil Engineering) for their insight and network layout data. 

Appendix A. FFactional-power scaled Chebyshev approximation. The 

numerical method frequently requires the height h as a function of the cross-sectional 
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Fig. 5.4. Case (I) results. Top: Initializing and optimal time series Q(0,t). Bottom: {dH/dx) 
over time for both inflow time series. 


area A. In the Preissman slot geometry, there is no analytic expression for h{A) 
below the pipe crown, necessitating some form of approximation. In what follows, all 
variables are assumed to be non-dimensionalized, e.g. A A/D^ and h —>■ h/D. 

One strategy is to use the relation A= g(0 — sin(0)) to solve for 9, and evaluate 
h = 5(1 — cos(0/2)). However, rootfinding is cumbersome as the computation must be 
performed for every cell during every time step. Furthermore, the authors found that 
representation of h{A) with a standard Chebyshev or other polynomial interpolant 
proved strikingly useless due to singular behavior at the origin. Comparing Taylor 
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Fig. 5.5. Case (II) results. Top: Initializing and optimal time series Q(0,t). Bottom: {dH/dx) 
over time for both inflow time series. 


series shows that h ^ near h = 0. Using a Chebyshev interpolant of the form 

N 

h{A) = Y,akTu(cA^/^-l), (A.l) 

k=0 

where Tk denotes the Chebyshev polynomial on [—1,1], scales out the singular 
behavior and ensures accuracy. The Chebyshev representation allows for fast eval¬ 
uation via a three-term recursion. The current version of the code uses iV = 20 
Chebyshev points for h{A), a change which speeds up the calculation of h{A) by a 
factor of 7 when compared to root-bracketing, without sacrificing double-precision 
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Table A.l 

Scaling parameter a values chosen. 


Function 

Range 

a. 

01 (A) 

0 < A < tt/S 

1/3 

02(A) 

7r/8 < A < At 

hjVl 

Ai(0) 

0 < 0 < 0 ( 71 / 8 ) 

1 

A2(0) 

0 ( 71 / 8 ) < 0 < 0(At) 

3/5 


accuracy. (Newton’s method is not desirable for this computation, as the derivative of 
f{9) = A — g(6* — sin(0)) is zero at the origin). 

We used the same strategy to better approximate the quantity 4>{A), which shows 
up in the speed estimation for the HLL solver and cannot be analytically expressed in 
terms of A. Several other authors use a ninth order McLaurin series for expressing (j) 
as a function of 9 [21, 16], which avoids repeated quadrature but still relies on root 
finding for 9 in terms of A. Furthermore, this expression also loses accuracy at the 
right end of the interval, with errors of order 10“^ near the pipe crown. 

The authors thus set out to improve both speed and accuracy in computation of 
(j) by finding a Chebyshev representation for both (f>{A) and = A{(j)), of the form 

N N 

cj){A) = aMcA^ - 1), A{c^) = ^ afcTfc(c</>“ - 1). (A.2) 

fe =0 fc =0 

As this geometry has peculiarities near both the pipe floor and pipe crown, we 
considered separate expansions on the first and second halves of the interval of interest. 
By numerical experiment with accuracy and coefficient decay of expansions, we found 
the best values of a for four expansions summarized in Table A.l. Figure A.l shows 
the accuracy and coefficient decay dependence on a for each of these cases. Our 
implementation uses N = 2Q terms in the expansions, but we tried several other 
values of N during our numerical experiments to make sure that that the number of 
Chebyshev points did not affect the scaling choices. 
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Fig. A.l. Coefficient decay and error as a function of power scaling a with N = 20 terms, 
for (piiA) (top) and Ai (bottom). The black hatched lines denote the expansion on the first half of 
interval (i = 1), and the green line denotes the expansion on the second half (i = 2). 
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